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Abstract 

It  has  long  been  established  that  supersonic  flow  over  axisymmetric  conical  bodies 
at  high  angles  of  attack  tend  to  develop  a  side  force  due  to  vortical  asymmetry.  One 
of  the  proposed  reasons  for  the  asymmetry  is  a  bifurcation  point  in  the  solution  of  the 
Navier-Stokes  equations.  This  study  investigated  the  possible  existence  of  a  bifurcation 
point  in  the  Navier-Stokes  equations  for  subsonic  laminar  flow.  Newton’s  method,  with 
gauss  elimination,  was  used  to  solve  the  steady-state,  viscous,  compressible  Navier-Stokes 
equations  in  spherical  cooordinates  assuming  conical  similarity. 
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A  NUMERICAL  DETERMINATION  OF  BIFURCATION  POINTS  FOR  LOW 


REYNOLDS  NUMBER  CONICAL  FLOWS 

I.  Introduction 

Asymmetric  vortices  on  the  leeward  side  of  slender  conical  bodies,  such  as  missiles 
or  fighter  aircraft  flying  at  high  angles  of  attack,  can  result  in  large  side  forces  on  the 
body,  even  at  zero  yaw.  These  large  side  forces  can  have  significance  in  vehicle  directional 
control.  The  development  of  symmetric  and  asymmetric  vortices  on  the  leeward  side  of 
slender  conical  bodies  has  been  observed  for  laminar  and  turbulent  flows,  subsonic  through 
hypersonic  speeds,  and  various  cross-sectional  shapes  (13). 


Experimental  studies  of  several  researchers,  including  (29),  (30),  and  (38),  have  identified 
four  distinct  flow  patterns  about  slender  bodies  at  angle  of  attack  and  zero-degree  sideslip. 
Using  6e  to  denote  the  cone  half-angle,  the  angle  of  attack  at  which  symmetric  vortices 
appear,  Oav  the  angle  of  attack  at  which  asymmetric  vortices  appear,  and  Out,  the  angle 
of  attack  at  which  unsteady  vort  'x  shedding  begins,  the  flow  about  a  slender  cone  can  be 
characterized  by: 

1)  At  low  angle  of  attack  (0  <  q  <  a,„),  axial  flow  dominates  and  the  flow  is  attached. 

2)  At  intermediate  angles  of  attack  (a,„  <  a  <  Oav),  a  symmetric  vortex  pair  appears 
on  the  leeward  side  of  the  body. 

3)  At  higher  angles  of  attack  (Oav  <  a  <  Oue),  crossflow  dominates  and  the  vortices 
become  asymmetric.  This  results  in  rapidly  increasing  side  forces,  even  at  zero  yaw. 

4)  At  very  high  angles  of  attack  (a„„  <  a  <  90®),  the  crossflow  dominates  completely, 
resulting  in  an  unsteady  vortex  shedding  similar  to  the  Karman  vortex  shedding  typical  of 
cylinders. 

Ericsson  and  Reding  (13),  and  Lowson  and  Ponton  (23),  report  that  for  laminar  flow, 
w  ^5,  Oav  «  20c,  and  au„  a  60®,  whereas  for  turbulent  flow  «  l.'iOc- 
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1.1  Historical  Background 

In  the  late  1940’s,  lateral  instability  caused  by  asymmetric  vortices  was  discovered 
in  flight  tests.  A  great  deal  of  experimental  work  followed,  including  Mudi  et  al.  (25), 
Peake  et  al.  (29),  Stahl  et  al.  (30),  Yanta  and  Wardlaw  (38),  and  Zilliac  et  al.  (39).  One 
of  the  original  computational  efforts  was  performed  by  Dyer  et  al.  (11),  who  showed  that 
asymmetric  vortices  on  the  leeward  side  of  a  slender  cone  could  be  predicted  by  standard 
line  vortex  models  using  symmetric  boundary  conditions.  This  study  was  followed  up  by 
a  more  detailed  study  by  Fiddes  and  Smith  (15),  and  Fiddes  (14),  where  vortex  sheets 
modeled  by  vortex  filaments  confirmed  the  earlier  results.  Assuming  conical  flow,  Marconi 
(24)  used  an  algorithm  based  on  the  Euler  equations  to  further  demonstrate  the  existence 
of  the  asymmetry  in  the  vortex  pair.  In  fact,  Marconi  found  that  at  a  critical  angle  of 
attack,  Oat) » the  symmetric  solution  was  unsteady.  He  observed  that  the  only  way  to  achieve 
a  symmetric  solution  was  to  impose  a  symmetric  condition. 

In  the  late  1980’s  and  early  1990’s,  more  powerful  computational  tools  brought  a 
virtual  explosion  to  the  area  of  calculating  asymmetric  vortical  flow.  Researchers  such 
as  Marconi  and  Siclari  (33),  Siclari  (32),  Degani  and  Levy  (9),  Degani  and  Schiff  (10), 
Degani  (7),  Kandil  et  al.  (19),  Batina  (3),  and  Vanden  and  Belk  (36)  are  but  a  few  of  the 
many  computational  efforts  in  this  area.  These  studies  range  from  subsonic  to  supersonic, 
laminar  and  turbulent.  Complementing  the  computational  efforts  were  many  experimental 
studies,  including  Degani  (8),  Zilliac  et  al.  (39),  Modi  et  al.  (25),  Yanta  and  Wardlaw 
(38),  Lowson  and  Ponton  (23),  Peake  et  al.  (29),  and  Stahl  et  al.  (34).  Many  of  these 
efforts  dealt  with  either  the  study  of  asymmetric  vortices  on  different  cross-sectional  shapes, 
including  delta  wings,  or  ihe  suppression  of  the  asymmetry  using  fins  or  strakes,  which 
impose  a  symmetry  in  the  flow. 

Several  of  the  numerical  solutions  utilized  the  conical  flow  assumption:  Batina  (3), 
Kandil  et  al.  (19),  Siclari  (32),  and  Siclari  and  Marconi  (33).  With  the  conical  flow  as¬ 
sumption,  changes  in  the  flow  variables  in  the  radial  direction  are  neglected,  so  that  the 
governing  equations  may  be  solved  in  two  dimensions.  The  advantage  of  using  this  assump¬ 
tion  is  that  a  much  more  detailed  numerical  analysis  can  be  performed.  The  disadvantage 
is  that  only  inviscid,  supersonic  flow  is  truly  conical,  where  the  length  scale  has  disap- 
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peared  from  the  problem.  Ericsson  (12)  investigated  experimental  results  for  supersonic 
and  subsonic  viscous  flow  and  concluded  that  "conical  flow  asymmetry  does  indeed  exist, 
but  only  up  to  moderate  angles  of  attack,  a  <  30  deg,  where  the  axial  flow  component  still 
has  a  strong  influence  on  the  crossflow-separation  characteristics.” 

All  of  the  results  from  the  studies  that  incorporated  the  conical  flow  assumption 
show  that  asymmetric  vortices  develop  around  the  critical  angle  Qo„  w  26^.  In  the  study 
performed  by  Siclari  and  Marconi  (33),  when  starting  from  a  symmetric  initial  condition 
the  residual  (difference  between  the  exact  and  approximate  solutions)  declined  by  ten 
orders  of  magnitude,  and  at  this  point  the  solution  was  essentially  symmetric  (Category 
2  above).  As  the  iterative  scheme  progressed,  the  residual  increased  to  almost  its  original 
value  and  then  declined  monotonically  to  machine  zero  where  asymmetric  vortices  were 
observed.  Monotonic  convergence  to  an  asymmetric  solution  could  be  achieved,  if  a  small 
asymmetry  was  introduced  into  the  initial  condition.  The  same  convergence  behavior  was 
also  seen  in  the  study  by  Kandil  et  al.  (19),  who  used  a  similar  algorithm.  These  studies 
utilized  a  time  integration  approach. 

Other  researchers,  (7),  (10),  and  (36),  computed  the  unrestricted,  three-dimensional 
flow,  over  a  slender  ogive-cylinder  body  at  angle  of  attack.  Vanden  and  Belk  (36)  computed 
both  supersonic  and  subsonic  flow.  In  both  cases,  a  localized  perturbation  to  the  body 
shape  was  needed  for  the  asymmetric  solution  (Category  3)  to  be  stable.  They  report  that 
this  observation  shows  that  vortex  *isymmetry  at  high  angles  of  attack  on  slender  bodies  is 
due  to  a  convective  instability  resulting  from  an  asymmetric  upstream  disturbance.  Degani 
(7)  computed  the  subsonic  laminar  flow  about  a  slender  ogive-cylinder  body.  Solutions  were 
obtained  for  angles  of  attack  ranging  from  a  =  20®  to  a  =  80®  and  a  Reynolds  number 
(based  on  freestream  conditions  and  cylinder  diameter)  oi  Re d  ~  200,000.  Results  at 
a  Mach  number,  M,  of  0.2  and  a  —  20®  showed  the  flow  to  be  steady  and  symmetric. 
The  introduction  of  a  space-invariant,  time-invariant  perturbation  placed  near  the  tip 
only  made  a  small  change  to  the  flow  structure.  At  a  =  40®  the  flow  was  steady  and 
symmetric,  but  became  asymmetric  with  the  introduction  of  the  perturbation.  The  level 
of  the  asymmetry  depended  on  the  size  and  location  of  the  perturbation.  When  the 
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perturbation  was  removed,  the  flow  returned  to  the  symmetric  case.  These  studies  also 
used  a  time  integration  approach. 

There  have  been  several  speculations  as  to  the  discrepancy  between  the  results  ob¬ 
tained  with  the  conical  flow  assumption  and  those  without  the  restriction.  In  regards  to 
the  results  obtained  by  Degani  and  Schiff  (10),  where  asymmetric  vortices  only  appeared 
when  a  perturbation  was  introduced  into  the  flow,  Kandil  et  al.  (19)  report  that  the  rea¬ 
son  for  this  is  the  result  of  ’’the  smallest  scale  of  the  grid  at  the  solid  boundary  and  the 
damping  effect  of  the  numerical  dissipation  in  the  axial  direction,  in  addition  to  the  grid- 
fineness  distribution.”  Vanden  and  Belk  (36)  claim  that  the  approximate  factorization  in 
the  numerical  scheme  used  by  Siclari  (32)  introduced  an  error  into  the  transient  solution. 

Complicating  the  comparitive  analysis  is  the  inability  to  achieve  a  perfectly  symmet¬ 
ric  test  in  experimental  work.  If  the  asymmetry  is  caused  by  arbitrarily  small  perturbations 
present  near  the  nose  tip  of  the  model,  then  the  asymmetry  will  be  observed  experimen¬ 
tally,  especially  since  machining  processes  used  to  construct  the  test  models  are  not  perfect, 
and  small  surface  imperfections  will  always  exist.  It  is  these  imperfections,  which  can  be 
large  compared  to  the  model  radius  near  the  model  tip,  that  could  be  responsible  for  flow 
asymmetry. 

1.1.1  Mach  number  effects.  Several  researchers,  including  (19),  (12),  (20),  and 
(13)  have  investigated  the  effect  of  Mach  number  on  the  angle  of  attack  at  which  the 
asymmetry  appears  and  the  relative  strength  of  the  vortices.  Their  conclusion  is  that  as 
Mach  number  increases,  ttav  increases,  while  is  not  affected.  This  trend  can  be  seen  in 
Figure  1.1,  from  (12).  Also,  the  relative  strength  of  the  side  forces  is  greater  at  subsonic 
speeds. 

1.1.2  Geometric  effects.  Again,  several  researchers,  (23),  (38),  (20),  (32),  and 
(13),  have  studied  the  effect  of  different  cross-sectional  shapes,  blunt  versus  sharp  noses, 
and  different  fineness  ratios  of  bodies.  The  cumulative  results  indicate  that  as  the  body 
becomes  thinner  (more  elliptic),  it  becomes  more  resistant  to  the  onset  of  the  asymmetry. 
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Figure  1.1  Mach  number  effects  on  Qou  (12) 
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Figure  L2  Possible  Solution  Spaces 

Siclari  found  that  asymmetric  vortices  develop  for  cross-sectional  shapes  other  than  circular 
cones.  His  study  included  elliptic,  biparabolic,  and  biwedge  cross-sections. 

Much  lower  side  forces  have  been  noted  for  blunt  versus  sharp  nosed  bodies.  This 
is  possibly  due  to  the  fact  that  a  bump  near  the  nose,  caused  by  surface  roughness,  has 
greater  impact  on  the  sharp  nose. 

i.S  Research  Objective 

The  main  objective  of  this  study  is  to  evaluate  the  solution  space  available  to  the 
Navier-Stokes  equations,  assuming  conical  similarity,  for  flows  around  a  slender,  circular 
cones  at  angle  of  attack.  A  direct  numerical  procedure  is  used  to  compute  bifurcation 
points  in  the  solution  space.  At  such  points,  symmetric  vortex  structures  become  unsta¬ 
ble,  while  stable  asymmetric  vortex  structures  become  admissable.  Most  of  the  previous 
investigations  in  this  area  have  dealt  with  the  supersonic  flow  regime.  This  study  provides 
additional  results  for  subsonic  flow. 
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Figure  1.2  shows  three  possible  solution  spaces  available  for  the  Navier-Stokes  equa¬ 
tions.  The  first  picture  (a)  represents  the  most  likely  solution  space,  and  the  one  presented 
by  Siclari  and  Marconi  (33).  In  these  figures,  the  ordinate  label  is  A,  which  represents  a 
free  parameter  in  the  solution.  For  this  analysis,  the  free  parameter  is  A  =  ajBc.  The  ab- 
sicca  represents  a  measureable  quantity  that  changes  after  the  bifurcation  point  is  passed. 
For  this  analysis,  the  side  force  coefficient  is  an  applicable  parameter.  The  trivial  path  is 
designated  as  the  path  where  the  side  force  coefficient  is  zero.  The  point  B  represents  the 
the  bifurcation  point.  To  the  left  of  B,  the  trivial  path  is  the  only  solution  available,  and  is 
considered  stable  (s).  To  the  right  of  B,  the  trivial  path  becomes  an  unstable  solution  (u), 
and  the  pitchfork  path  represents  the  stable  solution.  There  are  two  acceptable  solution 
paths,  one  above  the  ordinate  and  one  below.  These  solutions  are  identical,  mirror  images 
of  each  other.  For  this  analysis,  the  two  solutions  represent  which  side  of  the  cone  the 
asymmetry  occurs  on.  The  reason  Newton’s  method  with  gauss  elimination  was  chosen  as 
the  solution  technique  for  this  study  is  because  this  combination,  coupled  with  the  continu¬ 
ation  method,  allows  for  the  systematic  computation  of  the  entire  solution  space,  including 
the  unstable  branches.  Time  integration  techniques,  used  by  other  researchers,  can  not  be 
efficiently  coupled  with  a  continuation  method.  Also,  considering  the  convergence  behavior 
experienced  by  Siclari  and  Marconi  (33),  time  integration  routines  could  result  in  different 
solutions  of  the  computed  flowfield,  depending  on  the  convergence  criteria  used. 
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11.  Analysis 


In  this  Chapter,  a  model  problem  is  formulated  for  the  investigation  of  low-speed, 
conical  flows,  including  equations  of  motion,  boundary  conditions,  and  discretization.  The 
implications  and  validity  of  the  assumptions  of  the  model  formulation  are  covered  in  Sec¬ 
tion  2.1.  Sutherland’s  law,  used  to  model  the  coefficient  of  viscosity,  is  simplified  through 
linearization,  described  in  Section  2.2.  The  non-dimensional  governing  equations,  derived 
in  Appendix  A,  are  presented  in  Section  2.3,  while  boundary  conditions  are  described  in 
Section  2.4.  The  coordinate  transformation  from  physical  to  computational  space  is  out¬ 
lined  in  Section  2.5.  Spatial  discretization  of  the  physical  domain  is  covered  in  Section  2.6. 
The  approximation  of  the  derivative  terms  in  the  governing  equations  and  the  boundary 
conditions  is  shown  in  Section  2.7.  The  methodology  for  adding  fourth-order  numerical 
dissipation  is  presented  in  Section  2.8. 

S.l  Model  Assumptions 

The  assumptions  used  in  the  analysis  sire: 

•  Steady  flow 

•  Laminar  flow 

•  Conical  flow 

•  Thermally  and  caloricaJly  perfect  gas 

•  Constant  ratio  of  specific  heats 

•  Linearized  form  of  Sutherland’s  law 

•  Constant  Prandtl  number 

•  Stokes’  hypothesis 

•  Adiabatic  wall 

•  Normal  pressure  gradient  at  wall  is  zero 

Steady  flow  implies  |;^  =  0.  In  this  context,  and  throughout  the  remainder  of  this 
chapter,  /  represents  any  one  of  the  five  unknown  variables  in  the  problem:  density,  p,  the 
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three  velocity  components,  «,  v,  and  w,  and  internal  energy,  e.  Assuming  a  calorically  and 
thermally  perfect  gas,  implies 

P  =  pRT  €  =  c^T.  (2.1) 

The  ratio  of  specific  heats,  7,  is  assumed  to  be  1.4.  The  assumption  of  Stokes’  hypothesis 
means  that  the  bulk  viscosity  is  zero,  so  the  second  coefficient  of  viscosity.  A,  is  given  by 

A  =  -|/x. 

The  adiabatic  wall  assumption  implies  that  the  derivative  of  the  temperature  at  the 
wall  is  zero.  Linearizing  Sutherland’s  law  implies  a  small  variation  in  the  non-dimensional 
temperature  throughout  the  flowheld.  Figure  2.1  shows  the  comparison  of  Sutherland’s 
law  and  the  linearized  formula,  derived  in  Section  2.2.  Comparing  the  non-dimensional 
viscosity  given  by  Sutherland’s  law  (2.3)  and  the  linearized  formula  (2.10)  gives  a  difference 
of  2.99%  for  M  =  0.3  and  2.93%  for  M  =  0.7. 

With  conical  similarity,  changes  in  the  flow  variables  occuring  along  rays  (in  the  r 
direction)  emanating  from  the  cone  apex  are  assumed  to  vanish.  As  shown  in  Figure  2.2, 
r  is  one  of  the  three  coordinate  directions  for  a  spherical  coordinate  system.  Peake  and 
Tobak  (30)  illustrate  conical  flow  by  stating  that  "stream  surfaces  projected  on  to  concen¬ 
tric  spheres  centered  at  the  apex  (called  conical  flow  streamlines)  are  then  similar.”  The 
implications  are  that  all  derivatives  with  respect  to  r  of  the  flow  variables,  p,  u,  v,  w,  and 
e  vanish. 

Conical  similarity  is  exact  for  supersonic,  inviscid  flow  since  there  is  no  length  scale 
in  the  problem.  For  supersonic,  viscous  flow,  the  flow  is  not  exactly  conical,  because  in 
the  boundary  layer  the  flow  is  subsonic.  According  to  (30),  a  laminar  boundary  layer  is 
not  exactly  conical,  because  the  boundary  layer  grows  as  r°  and  in  a  fully  turbulent 
boundary  layer,  where  the  boundary  layer  grows  as  r,  the  flow  is  nearly  conical.  Peake 
and  Tobak  (30)  summarize  the  issue  of  conical  flow  by  stating: 

When  the  Reynolds  number  is  sufficiently  high  so  that  transistion  occurs  in 
proximity  to  the  apex,  the  near-conical  nature  of  the  experimentally  measured 
flow  demonstrates  a  virtual  absence  of  length  effects  in  the  streamwise  direc¬ 
tion:  the  flow  is  dominated  completely  by  the  circumferential  pressure  field. 
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Thus,  the  characteristics  of  these  flow  flelds  can  be  determined  through  mea¬ 
surement  or  by  computation  at  essentially  one  streamwise  station.  In  fully 
turbulent  and  fully  laiminar  subsonic  freestream  flow,  even  though  base  and 
thickness  effects  become  measureable,  the  circumferential  pressure  gradients 
stiU  dominate,  to  the  extent  that  virtual  conicity  of  the  separation  lines  and 
shear-stress  directions  are  still  maintained. 

The  applicability  of  the  conical  flow  assumption  for  subsonic  flows  is  the  subject  of  a 
paper  by  Ericsson  (12).  In  this  paper  Ericsson  reviewed  experimental  results  and  concluded 
that  "conical  flow  asymmetry  does  exist  on  very  slender  cones  because  of  the  still-present 
strong  axial  flow  component  at  26c  <  a  <  30*.”  On  cones  where  the  asymmetry  does  not 
develop  until  a  >  30*,  the  flow  is  nonconical.  For  this  reason,  the  cone  half-angle  used  in 
this  study  is  =  5  degrees. 

The  assumptions  of  a  thermally  and  calorically  perfect  gas,  along  with  constant 
ratio  of  specific  heats,  is  valid  for  the  subsonic  and  low  supersonic  flows  investigated  in 
this  study  (2).  Reference  (37)  supports  the  assumption  of  constant  Pr  and  zero  pressure 
gradient  normal  to  the  wall  in  the  boundary  layer.  Symmetric  and  asymmetric  vortices 
have  been  shown  to  exist  in  flows  over  circular  cones  for  steady,  laminar  flow  (References 
(13),  (19),  (33),  (7)).  No  attempt  was  made  to  model  turbulent  flows. 

2.2  Linearization  of  Sutherland’s  Law 

Since  viscosity  is  not  assumed  to  be  a  constant,  a  relationship  is  desired  that  expresses 
viscosity  in  terms  of  one  or  more  of  the  unknown  variables.  Sutherland’s  law  (1)  provides 
this  relationship: 

m(T)  =  (2.2) 

where  Cj  =  llO.d'iif  and  C\  will  be  eliminated  through  non-dimensionalization.  Using 
(2.1),  fi  can  be  expressed  as  a  function  of  the  internal  energy,  e.  However,  with  this 
expression,  the  complexity  of  the  derivation  of  analytical  Jacobian  elements  is  greatly 
increased,  (examples  of  the  derivation  of  Jacobian  elements  appear  in  Appendix  B).  To 
avoid  such  complexity,  (2.2)  is  first  placed  in  non-dimensional  form  and  then  linearized  to 
provide  a  simpler  expression  for  fi. 
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2.00 


Temperature 

Figure  2.1  Comparison  of  Sutherland’s  law  and  linearized  formula 


The  freestream  viscosity  is  found  by  evaluating  (2.2)  at  the  fireestream  temperature. 
Too-  The  appropriate  non-dimensional  scale  factors  for  n  and  T  are  then  the  freestream 
values,  /Xqo  ^d  Too  (sU  of  the  non-dimensional  scale  factors  are  given  in  Appendix  A). 
Using  the  superscript  (*)  to  denote  a  non-dimensional  value,  equation  (2.2)  is  written  in 
non-dimensional  form  as 

*  _  Jfill  _  (C2  +  T00)  _  C2/T„  +  1 

^  (C2-FT)  \T^J  C2lT^+TIT^- 

Using  T*  =  T/Too  and  cj  =  Cj/Too,  (2.3)  is  rewritten  as 

=  (rf'"  .  (2.4) 

Equation  (2.4)  is  linearized  by  expanding  it  in  a  first-order  Taylor  series  expansion  about 

=  m*(T4)  +  AT*  -f  0(AT*)T  (2.5) 
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where  e  and  T  are  non-dimensionalized  using  the  scale  factors  provided  in  Appendix  A, 
resulting  in 

RT'*T' 

=  =>  r*  =  7(7-l)MV.  (2.9) 

Substituting  (2.9)  into  (2.8)  results  in  a  linear,  non-dimensional  form  of  n  in  terms  of  the 
internal  energy: 


H*  =  cie*  -I-  Ca, 


where  Ci  and  Ca  are  given  by: 


Cl 


1  ' 
cj-f  1. 


7(7-  1)M^ 


1  • 
c;-!- 1. 


(2.10) 


(2.11) 


The  thermal  conductivity,  k,  is  related  to  through  the  assumption  of  a  constant  Prandtl 
number,  Pr,  where 


K 

Non-dimensionalization  of  (2.12)  yields 


Pr  = 


Cp/X>oo 

K*Koo 


When  Pt  =  constant,  (2.13)  provides 


(2.12) 


(2.13) 


H*=k\  (2.14) 

Equation  (2.10)  is  used  to  replace  the  non-dimensional  thermal  conductivity,  k*,  which 
appears  in  the  energy  equation  (derived  in  Appendix  A)  For  the  remainder  of  the  analysis, 
the  superscript  (*),  representing  a  non-dimensional  variable,  is  dropped  for  convenience. 


2.3  Governing  Equations 

Following  the  assumptions  outlined  in  Sections  2.1  and  2.2,  the  non-dimensional 
governing  equations  are  developed  in  Appendix  A.  A  schematic  of  the  spherical  coordinate 
system  is  given  in  Figure  2.2.  With  the  subscripts  9  and  <j)  denoting  differentiation,  the 
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non-dimensional  governing  equations  in  spherical  coordinates  are: 
Continuity  equation 


2pu  sinO  pv  COS0  -I-  vpg  sin^  -f  pvg  sin0  +  -|-  wpg,  =  0 

r-momentum  equation 


pr 


.  2  2 
+  ~r-^  —  V  —  w 

sinO 


,  r  1  ,  CiC^ 


w 

sinff  sin^5 


Ci€ -|- C2  [SU  7uCOtd  IVg  7tVg, 


Re 


'8u 

.T 


-  “««  -  cot^ 


3  sind  sin^d 


=  0 


0-momentum  equation 


pr 


vvg  -I-  +UV  -w^  cotdj  -f-  r(7  -  1)  [pe#  -|-  ep,] 


c\ee 

'2v  coi0 

2w^ 

Ave 

2u 

'wco&0 

We 

Vg, 

Re 

3 

3sind 

~  T 

~  T. 

Re 

.  sin^^ 

sind 

sin^d. 

cie  -I-  C2 

Re 


3  sin  5 


wgg,  8ue  4vj  cot^  2v 

3sin^  3  3  3sin^5 


3  sm^0 
=  0 


4>-momentum  equation 


pr 


WWg, 

vwe  -i — r— T  -I-  wu  +  vw  cot0 
sin0 


Cie« 


Re 

Cie  +  C2 


We  —  w  cot0  -I- 

smS 


Cie^ 

Re 


2u  4v  cosO 

+  „  .  X.  + 


2ve 


3  sin^C  3  sin0  3  sin^^  3  sind 


Re  l3  sin^^ 


,  ...  ,  '^e4>  .  8w^  ,  Tv^cos^ 


3  sintf  3  sin^  3  sin^^ 


w 


2w  "j-  -I-  We  cot0  -  2w  cot^0 
sm  0 


=  0 


(2.15) 


(2.16) 


(2.17) 


(2.18) 


Energy  equation 


7(cie  -f-  C2) 

PrRe 


eg  cot9  +  e*#  +  — 


ICi 


sin^dj  PrRe 


^1  + 


sin^O 


CiC  +  C2 

Re 


4vg  4uvg  4u^  4wl  4t;^cot^^ 
— Sl  + 1  + J. 2_  4. - 

3  3  3  SsinH  3 


4uwj,  Sutu^cotd  4ut;cotd  v?  ,9^ 

+ - ^  + - - - + - +  — tr-  i-wl  +  cot^e 

3sin^  3sind  3  sin^^  * 

,  2wgVs  2wv^cot6  „  ul  2wu^ 

sm(9  Sind  sin^d  sind 


4-w^  +  +  uj  —  2vug  — 


4vgw^  4vvg  cotdl 


3sind 


=  0 


(2.19) 


2.4  Boundary  Conditions 

The  governing  equations  are  solved  at  all  the  internal  nodes,  therefore  only  two 
physical  surfaces  require  boundary  conditions.  These  surfaces  are  shown  in  Figure  2.2, 
where  the  inner  circle  represents  the  wall  boundary,  and  the  outer  circle  represents  the 
freestream  boundary.  For  the  wall,  the  pressure  derivative  normal  to  the  wall  is  set  equal 
to  zero,  and  an  adiabatic  wall  condition  is  used  for  aU  the  runs.  This  implies. 


(2.20) 


(2.21) 


These  boundary  conditions  are  used  in  all  the  other  computational  efforts  in  this  area  (7), 
(10),  (19),  (32),  (33),  and  (36).  The  normal  pressure  derivative  can  be  considered  a  first 
step  approximation,  where  a  more  accurate  boundary  condition  would  be  setting  normal 
momentum  equal  to  zero.  Using  equation  (A. 12),  which  relates  pressure  to  density  and 
internal  energy  using  the  perfect  gas  assumptions,  along  with  equation  (2.9),  which  relates 
temperature  to  internal  energy,  equations  (2.20  -  2.21)  can  be  expressed  as 


1  =  0 

(2.22) 

de 

(2.23) 
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Wall 

Freestream  Boundary 

P=  1 

U  =  Uf, 

u  =  0 

V  =  Vj, 

111  =  0 

w  =  0 

e  =  ej. 

Table  2.1  Boundary  conditions  for  supersonic  flow 

For  the  viscous  flow  cases,  no  slip  at  the  wall  is  enforced;  u  =  v  —  w  =  0.  For  the  inviscid 
flow  cases,  slip  flow  and  impermeability  give  |^  =  Oj  |^  =  0,  and  t;  =  0. 

On  the  freestream  boundary  (denoted  by  the  subscript  fs),  a  general  set  of  equations 
is  used  to  determine  the  values  for  the  flow  velocities  in  terms  of  the  angle  of  attack,  a, 
and  the  coordinate  angles,  8  and  0: 

«i  =  cosacosfl  U2  =  sin  a  sin  d  cos  ^  uj,  =  {u\  +  (2.24) 

I’l  =  -  cos  a  sin  0  i;2  =  sin  a  cos  cos  vj,  =  {v^  +  (2.25) 


tni  =  0  ti;2  =  -  sin  a  sin  </)  wj,  =  {wl (2.26) 

Equations  (2.24  -  2.26)  are  also  used  to  initialize  the  flow  velocities. 

An  expression  for  non-dimensional  internal  energy  at  the  freestream  boundary,  ej, , 
is  obtained  by  evaluating  equation  (2.9)  at  Tj,  =  1.  Therefore  (dropping  the  (*)  for 
convenience). 


1 

The  boundary  conditions  for  each  case  are  summarized  in  Tables  (2.1  -  2.3).  For  sub¬ 
sonic  flow,  the  continuity  equation  is  used  as  a  freestream  boundary  condition.  Section  2.6 
outlines  the  use  of  the  continuity  as  a  boundary  condition. 


2. 5  Coordinate  Transformation 

In  computational  fluid  dynamics,  the  governing  equations  and  boundary  conditions 
are  solved  at  a  finite  number  of  discrete  points  which  represent  the  solution  domain.  In 
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Table  2.2  Boundary  conditions  for  subsonic,  invisdd  flow 


WaU 

Freestream  Boundary 

continuity 

=  1 

V  =  0 

P=  1 

w  =  Wf, 

e  =  Cy, 

Table  2.3  Boundary  conditions  for  subsonic,  viscous  flow 


Figure  2.3  Physical  and  Computational  Domains 
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the  physical  (</>,  ff)  space,  grid  points  need  to  be  clustered  near  the  wall  so  the  boundary 
layer  can  be  properly  resolved.  This  can  be  achieved  by  either  uneven  spacing  between 
grid  points  in  the  radial  direction,  or  by  placing  many  mora  evenly  spaced  points  along 
each  radial  line.  Clearly  the  first  approach  is  more  efficient  computationally.  However, 
the  standard  second-order-accurate  difference  approximations  used  in  the  discretization 
of  the  governing  equations  require  a  uniform  spacing  between  nodes.  Even  spacing  is 
obtained,  with  node  clustering  in  the  physical  domain,  through  a  coordinate  transformation 
of  the  physical  domain  to  a  computational  domain.  In  this  problem,  the  physical  space 
is  designated  by  the  coordinate  pair  and  the  computational  space  is  designated  by 

the  coordinate  pair  (^,q),  where  Figure  2.3  shows  the  mapping  between  the  two  domains. 
For  this  study,  the  physical  domain  could  be  separated  anywhere,  but  the  windward  side 
of  the  cone  was  selected  because  the  gradients  in  the  dependent  variables  are  smaller  in 
this  region.  Using  subscripts  to  denote  differentiation,  the  transformation,  from  (16),  is 
expressed  in  general  terms  as 

=  e  =  e{i,Ti),  (2.28) 

which,  from  the  chain  rule  of  partial  differentiation,  leads  to 

d<j)  =  d6  =  6(d^  -f-  9,,dri.  (2.29) 

Similarly,  the  reverse  transformation  yields 

^  =  r,  =  q(0,6»),  (2.30) 

d^  =  i^d4>  -I-  ^edO,  dr]  =  T]^d4>  +  rjgdd.  (2.31) 

When  expressed  in  matrix  form,  the  transformation  is 


6 

1 _ 

<i>ri 

-1 

1 

0^ 

Vi 

. 

“  7 

-0, 

4>(. 

where  the  Jacobian  is  defined  as 


J  — 


2-11 


The  metrics  of  the  transformation  are 

c  -  h.  c  -  Ti  -  n  =  ^ 

U  j  j  j  j- 

The  metrics  and  the  Jacobian  are  evaluated  numerically  at  all  node  points  in  the  domain 
using  the  difference  approximations  provided  in  Section  2.7. 

The  spatial  derivatives  are  transformed  from  physical  to  computational  space  with 
the  following  formulas,  from  (16)  and  (35), 

f<t>  =  +  ’?«/»)>  ft  =  6/?  +  Veftii  (2.32) 

fu  =  (Ifu  +  +  vlfnv 

-  +  vl<t>r,r,){Uh  +  (2-33) 

fet  =  Chi  +  26t)tf/ei,  + 

“  +  2C«»7»0{i,  +  V9^rm)i^tf(  +  Vtfft) 

~  +  ‘n]4>nr)){i<t>h  +  (2.34) 

f9<t>  =  +  {Vt^<t>  +  it‘n<t>)hn  +  V^Vsft},, 

+  {Ve^vv  +  ^t^iv  ~ 

~  +  Ve^ifi  +  VeV^i’fn  +  (2.35) 

~j  =  +  Ve^(ri  +  'n<i><(>(ri  +  (2.36) 

-J  =  +  Ve^rtn  +  'n<t><t>vr)  +  ^t^(v  (2.37) 
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Figure  2.4  Nine-point  Stencil  for  interior  nodes 

2.6  Discretization  of  the  Governing  Equations  and  Boundary  Condi*  is 

Following  the  transformation  of  the  governing  equations  to  the  (^,  q)  computational 
space,  the  spatial  derivatives  are  approximated  to  second-order-accuracy  with  discrete 
central-difference  formulas.  First  and  second  derivatives  of  ^  and  q,  as  well  as  cross¬ 
derivative  terms,  appear  in  the  transformed  equations.  Application  of  central-difference 
formulas  at  a  given  node,  k,  requires  information  from  adjacent  node  points.  For  the 
central-difference  formulas  required  in  this  study,  the  influential  neighbors  of  an  interior 
node  point  are  contained  within  a  nine-point  computational  stencil,  shown  in  Figure  2.4. 
The  difference  approximations,  from  (1),  used  for  the  interior  nodes,  are 


r  1  fke  fkw 

'  2AC  ’ 

(2.38) 

t  1  fkn  ~  fks 

“  2Av  ’ 

(2.39) 

el  fke  -  ^fk  +  fi-j, 

m\k  - 

(2.40) 

,  1  fkn  -  2/*  +  fk$ 

Mk  - 

(2.41) 

el  1  fkne  “  fk$e  fknw  fksw  \ 

~  2A^  V  2Arj  2At}  )  ' 

(2.42) 
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Figure  2.5  Stencil  for  one-sided  boundary  derivatives 


For  the  wall  boundary  nodes,  normal  derivatives  are  replaced  directly  with  ij  deriva¬ 
tive  terms.  One-sided,  second-order-accurate  difference  approximations  are  used  to  ap¬ 
proximate  the  T)  derivatives: 

f  —  ~^/t  •h  ^hn  —  fknn 

~  2A| 

where  Figure  2.5  shows  the  stencil  used  for  the  one-sided  differences  needed  for  the  bound- 
aury  derivatives.  A  square  node  arrangement  is  used  for  all  cases,  so  =  A?7  =  1. 


2. 7  Implementation  of  Continuity  at  the  Freestream  Boundary 

For  subsonic  flow,  the  continuity  equation  (2.15)  is  used  as  a  freestream  boundary 
condition.  This  equation  is  written  as 

2pusin^  -t-  (pt7sin^)j  -f  {pw)^,  =  0  (2.44) 

The  equation  is  rewritten  in  this  form  so  that  the  first  and  second  terms  can  be 
evaluated  at  the  midpoint  between  the  boundary  node  and  the  next  node  in  from  the 
boundary,  shown  in  Figure  2.6.  The  difference  approximation  (2.36)  is  used  to  represent 
the  ^  derivative  term.  The  first  term  is  represented  by  averaging  the  values  at  the  midpoint. 


Figure  2.6  Stencil  for  Continuity  equation  at  freestream  boundary 


A  second-order-accurate  difference  about  the  midpoint  is  used  to  represent  the  r]  derivative 
term,  so  only  two  node  points  are  required.  The  approximation  equation  is 


fk  -  fk, 

Ati 


(2.45) 


All  of  the  grids  used  in  the  study  have  evenly  spaced  radial  lines  around  the  cone,  covered 
more  extensively  in  section  2.7  .  This  results  in  the  metrics,  ^0  and  being  zero.  Therefore 
the  continuity  equation  used  for  the  freestream  boundary  equation  is 


(pusinfl)*  -I-  (pusintf)*,  •+•  i]0  [(pusind)*  -  (pvsind)*,] 

+  Y  ”  (/»"')*«-]  =  0*  (2-46) 

2.8  Grid  Generation  and  Discretization 

Two  different  grids  are  used  in  the  analysis.  The  first  is  used  for  the  inviscid  cases, 
both  supersonic  and  subsonic.  For  these  cases,  there  is  no  boundary  layer  to  resolve  so 
there  is  no  need  to  cluster  the  points  near  the  cone  surface.  Therefore,  the  radial  lines 
are  evenly  spaced  around  the  cone,  and  the  points  are  evenly  spaced  along  each  radial 
line.  When  the  flow  is  viscous  and  subsonic,  node  points  must  be  clustered  near  the  cone 
surface  to  accurately  resolve  the  boundary  layer.  The  radial  lines  are  still  spaced  evenly 
around  the  cone,  but  the  nodes  along  each  radial  line  are  distributed  using  a  geometric 
progression,  outlined  in  (6).  Figures  (2.9)  and  (2.10)  show  examples  of  the  different  grids. 

The  nodes  in  the  computational  domain  are  numbered  using  a  row-by-row  convention, 
as  shown  in  Figure  2.7.  The  advantage  of  numbering  the  nodes  this  way,  as  compared  to 
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Node 

II 

i  =  I 

2<t</-l 

ke 

ifc+1 

k-I+1 

*  +  1 

kw 

k  +  I-l 

k-l 

ifc-1 

kn 

k  +  I 

k  +  I 

k  +  I 

ks 

k-I 

k-I 

k-I 

kne 

k  +  I+l 

k  +  1 

k  +  I  +  l 

knw 

k  +  2I-l 

k  +  I-1 

k  +  I-1 

kse 

k-I+l 

k-2I+l 

k-I+1 

ksw 

Ar  —  1 

_ 1 

k-I+l 

k-I-1 

Table  2.4  Node  relationships  for  nine-point  stencil 


column-by-column,  is  that  the  governing  equations  could  be  used  with  the  standard  nine- 
point  stencil,  shown  in  Figure  2.4,  at  the  nodes  adjacent  to  the  separation  region,  shown 
in  Figure  (2.3).  If  the  nodes  are  numbered  in  a  column- by-column  scheme,  then  boundary 
equations  would  be  necessary  for  the  nodes  bordering  the  separation  region.  An  example  of 
this  is  the  two-dimensional  flow  over  a  circular  cylinder,  where  the  grid  is  typically  cut  on 
the  windward  (upstream)  side  of  the  cylinder  and  symmetry  boundary  conditions  are  used 
for  the  node  points  at  the  cut.  This  situation  is  avoided  in  this  study,  because  symmetry 
boundary  conditions  may  affect  the  development  of  asymmetric  vortices.  The  disadvantage 
of  row-by-row  numbering  is  the  bandwidth  of  the  Jacobian  matrix  is  directly  related  to  the 
second-order-accurate  difference  approximations  used  to  approximate  the  wzJl  derivative 
boundary  conditions,  and  cross-derivative  terms  in  the  momentum  equations. 

Letting  i  represent  the  radial  line,  and  j  represent  the  node  on  a  radial  line,  starting 
with  1  at  the  wall,  the  node  number  k  is  found  by 


k  =  (j-l)*I  +  i.  (2.47) 

Table  2.8  describes  the  relationship  between  the  primary  node  k  and  the  nodes  in 
the  nine-point  stencil  at  different  radial  lines  in  the  domain. 

As  can  be  seen  from  Table  2.8,  the  bandwidth  caused  by  numbering  the  nodes  using 
a  row- by-row  scheme  is  4I-I-1  due  to  using  second-order  accurate  one-sided  differences  to 
represent  the  wall  boundary  conditions.  The  maximum  bandwidth  needed  for  the  sepa¬ 
ration  nodes  (i  =  1,  i  =  /)  is  41-1.  The  computational  speed  for  the  solution  algorithm, 
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Figure  2.7  Node  numbering  in  computational  domain 

outlined  in  Chapter  III,  is  largely  determined  by  the  bandwidth  of  the  Jacobian  matrix. 
As  the  bandwidth  is  increased,  the  time  required  for  convergence  increases  exponentially. 

2,9  Added  Numerical  Dissipation 

For  the  subsonic,  viscous  flow  cases,  fourth-order  numerical  dissipation  is  added  to 
provide  numerical  smoothing.  This  smoothing  diminishes  the  very  small  oscillations  in 
the  computed  flow  variables  present  in  the  solutions  near  the  freestream  boundary,  caused 
by  a  lack  of  grid  resolution.  The  added  dissipation  also  aflects  the  numerical  results 
near  the  body,  so  care  must  be  taken  when  determining  the  amount  of  extra  dissipation 
added.  Figure  2.8  shows  the  stencil  needed  for  the  added  dissipation.  The  difference 
approximations,  from  (1),  for  fourth-order  spatial  derivatives  are 

~  ifkee  -  ^fke  +  6/*  “  4/jfc«,  -f-  fkww)  ,  (2.48) 

*  (Ann  -  4fkn  +  ^fk  “  ^fks  +  A**)  >  (2.49) 

where  is  an  adjustable  parameter. 
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60x81  Packed  Grid 


Fi»uro  2.10  Packed  Grid 
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III.  Solution  Algorithm 


3.1  Newton’s  Method 


Newton’s  method  is  an  iterative  scheme  that  solves  the  nonlinear  system  of  equations  of 
the  form 


Fix)  =  0,  (3.1) 

where  *  is  a  vector  of  N  unknowns,  and  F  is  a  set  of  N  equations.  The  Newton  iteration 
formula  is 


F;(r+i-r)  =  -F(x‘'), 


where  is  the  Jacobian  matrix,  whose  elements  are  given  by 


(3.2) 


SF, 

dF, 

SF, 

9*1 

0*3 

a*. 

a£i 

OFa 

Ml 

dxi 

0*3 

9x, 

OF. 

OF. 

.  .  Mm. 

dxi 

0*3 

9*« 

(3.3) 


x"  represents  a  known  solution  state  that  is  improved  with  the  solution  of  (3.2),  Ax  = 
x""*"^  -  x".  Using  central  differences  and  a  nine-point  stencil  to  represent  the  derivatives 
in  the  governing  equations,  (3.2)  results  in  a  banded  system  that  is  solved  using  Gauss 
elimination.  The  solution  of  the  system  is  Ax  =  x""*"*  —  x".  After  each  iteration,  the 
approximate  solution  is  updated  by  i*"''*  =  Ax  -I- x".  Successive  iterations  are  computed 
until  the  largest  absolute  Ax  is  less  than  some  small  value  c,  or  reaches  machine  zero. 
Newton’s  method  is  guaranteed  to  converge  quadratically  if  the  Jacobian  matrix  is  non¬ 
singular  and  the  initial  guess  is  sufficiently  close  to  the  exact  solution  (18). 


The  analytical  Jacobian  elements  are  determined  by  differentiating  each  of  the  equa¬ 
tions  with  respect  to  the  unknown  variables.  An  example  of  how  these  elements  are  derived 
is  shown  in  Appendix  B.  For  the  discrete  equations  in  this  study,  the  Jacobian  elements  are 
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relatively  complicated.  There  is  always  the  chance  that  there  are  errors  in  the  derivation 
and  implementation  of  the  elements.  There  are  two  methods  of  checking  the  accuracy  of 
these  elements.  The  first  is  if  the  scheme  is  quadratically  convergent,  then  the  Jacobian 
matrix  is  considered  to  be  correct.  The  second  method  is  to  compare  the  analytical  Ja¬ 
cobian  elements  with  elements  developed  numerically.  Using  a  method  implemented  by 
Morton  (27),  the  elements  of  the  two  matrices  are  compared  numerically,  and  when  they 
are  equal  the  analytical  elements  are  considered  correct.  This  numerical  scheme  consists 
of  solving  the  function  evaluation  (right-hand  side)  for  a  given  solution  vector,  x,  at  each 
node,  and  then  perturbing  the  solution  vector  by  a  small  amount,  i  +  f ,  and  then  solving 
the  function  evaluation  with  the  perturbed  solution  vector.  Subtracting  the  two  results 
for  each  unknown  at  each  node  in  the  nine-point  stencil  results  in  the  numerical  Jacobian 
elements.  Since  a  finite  perturbation  is  used  (c  =  0.000001),  the  accuracy  of  the  numerical 
elements  is  taken  to  be  one  order  less  than  the  t  used.  Use  of  this  method  allows  for  a 
systematic  way  of  checking  each  of  the  Jacobian  elements,  since  any  differences  in  the  two 
Jacobian  matrices  leads  directly  to  the  probable  incorrect  element.  This  allows  the  analyt¬ 
ical  Jacobian  matrix  to  be  verified  with  confidence  in  a  short  period  of  time.  Without  this 
type  of  comparison,  there  is  always  the  possibility  that  incorrect  development  or  imple¬ 
mentation  of  the  analytical  Jacobian  elements  is  the  cause  of  any  difficulties  in  obtaining 
a  correct  solution. 

3.2  Continuation  Method 

It  has  been  proven  that  if  a  solution  x*  is  known  and  is  nonsingular,  then  for 
some  range  of  A  about  A*  there  exists  a  unique  solution  path  through  (x*,A*).  The  proof 
is  outlined  in  (21).  Pseudo- arclength  continuation  (PAC)  is  Keller’s  method  of  computing 
solutions  along  the  solution  path  by  using  information  at  (x*,A*).  to  compute  the  next 
solution  point. 

Figure  3.1  is  representative  of  a  solution  path  found  by  plotting  the  norm  of  the 
solution  vector  versus  the  free  parameter  A.  The  PAC  process  is  to  compute  a  tangent 
vector  r  at  a  known  solution  i* ,  designated  by  P.  Then  at  a  distance  d  away,  search  along 
a  line  perpendicular  to  T  for  the  next  solution  point.  This  is  done  by  first  using  arclength 
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to  parameterize  the  solution  path  {x  —  x(s),X  =  A(5),  and  F  =  F{s)  =  0).  Then  the 
tangent  vector  is  computed  by 

^F(x(s);A(s))  =  0.  (3.4) 

Using  the  chain  rule,  equation  3.5  becomes 

F^{s)ic{s)  +  F,{s)X{s)  =  0,  (3.5) 

where 

j(s)  =  g(s),  (3.6) 

and 

A(s)  =  £(i).  (3.7) 

The  definition  of  arclength  is 


llil|2  +  A^(s)  =  1.  (3.8) 

Equations  3.5  and  3.8  can  then  be  solved  as  a  system  for  the  tangent  vector: 


providing  the  Jacobian  matrix  is  not  singular.  Now  define  4>  such  that 


ci>  =  f;'f,{s). 

Then  equations  3.5  and  3.8  give  the  following  relationships 


A(s) 

and 


±1 

yiTpiF’ 


(3.10) 


(3.11) 


a:(s)  =  —X{s)<f>.  (3.12) 

The  sign  of  equation  3.11  represents  the  direction  of  the  tangent  vector  and  is  therefore 
indeterminate.  At  the  startup  of  the  continuation  process  this  sign  is  set  depending  on 
which  part  of  the  solution  path  is  to  be  computed. 
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FVom  the  solution  P,  the  tangent  vector  T,  and  the  distance  d,  the  initial  solution 
vector  Qo  be  computed: 


The  solution  Q  lies  on  a  line  perpendicular  to  T  passing  through  Qq.  This  condition  can 


be  stated  mathematically  as 


D  =  ipilSQ-Jp)  +  (XQ- Xp)Xp  =  d  (3.14) 

and  can  be  added  to  the  system  of  nonlinear  equations  3.1.  This  new  system  is  then  solved 
by  Newton’s  method  with  Qo  as  the  initial  guess.  Qo  becomes  a  better  and  better  first 
approximation  as  d  gets  smaller.  The  solution  Xq  is  obtained  when  F  =  0  and  D  =  d. 
Another  r  '  ion  can  be  computed  by  repeating  the  process. 

The  entire  section  on  the  continuation  method  was  taken  directly  from  Morton  (26). 
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IV.  Navier-Stokes  Solutions 


Initially,  the  boundary  conditions  in  Table  4.1  were  used  to  solve  the  subsonic  flow 
over  a  slender  cone  at  zero  angle  of  attack.  These  boundary  conditions  are  the  same  as 
those  used  by  Degani  and  SchifF  (10)  in  their  three-dimensional  analysis.  However,  when 
these  boundary  conditions  are  used,  oscillations  (odd-even  decoupling)  developed  in  the 
computed  flow  variables  p,  u,  u,  w,  and  e.  A  significant  amount  of  effort  was  spent  trying 
to  find  an  error  in  the  code  that  would  account  for  an  error  of  this  type.  The  derivation  of 
the  governing  equ  -ions  was  checked  by  two  independent  methods,  outlined  in  Appendix 
A.  The  derivation  of  the  analytical  Jacobian  elements  and  their  implementation  in  the 
code  was  checked  using  the  numerical  procedure  outlined  in  Chapter  III.  The  placement 
of  the  Jacobian  elements  into  the  matrix  was  checked  extensively,  especially  in  the  grid 
separation  region.  Finally,  a  comparison  was  made  with  a  1-D  inviscid  code  at  o  =  0,  using 
the  boundary  conditions  in  Table  4.2.  These  boundary  conditions  represent  slip  flow  and 
impermeability  at  the  wall  boundary.  The  two  methods  gave  the  same  numerical  result, 
both  having  oscillations  in  the  computed  solution.  The  addition  of  fourth-order  numerical 
dissipation,  outlined  in  Chapter  2.9,  did  not  reduce  the  oscillations.  Next,  the  boundary 
conditions  in  Table  4.2  were  used  for  a  supersonic  inviscid  model  at  a  =  0,  and  the  solution 
matched  that  of  exact  methods,  both  in  shock  angle  and  velocity  components.  The  results 
of  this  analysis  is  presented  in  Section  4.1.  This  analysis  helped  to  verify  the  correctness 
of  some  of  the  computer  code,  especially  the  implementation  of  the  Jacobian  elements  in 
the  separation  region.  At  this  stage,  the  freestream  boundary  conditions  were  suspected 
as  the  cause  of  the  oscillatory  behavior  in  the  solution. 

Several  attempts  were  made  to  develop  freestream  boundary  conditions  for  which 
the  solution  would  be  free  of  oscillations.  Several  variations  were  attempted,  and  the 
boundary  conditions  outlined  in  Table  4.3  and  Table  4.4  were  found  to  give  good  results. 
The  key  to  the  problem  is  how  the  continuity  equation  is  implemented  at  the  discrete  node 
points.  Initially,  the  continuity  equation  was  used  as  a  freestream  boundary  condition,  as 
well  as  one  of  the  governing  equations  for  the  interior  nodes.  This  method  still  resulted  in 
oscillatory  behavior.  Then,  the  continuity  equation  was  evaluated  at  the  midpoint  between 
the  node  in  question,  k,  and  the  next  node  closer  to  the  wall,  ks  (see  Figure  2.6),  for  both 
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boundary  and  interior  continuity  equations.  This  allows  the  use  of  second-order-accurate 
difference  approximations  for  the  freestream  boundairy  derivatives  while  only  using  two 
nodes  instead  of  three.  Good  results  were  obtained  depending  on  the  number  of  grid 
points  along  each  radial  line  relative  to  the  difference  between  the  cone  haJfangle,  0ci  ^md 
the  angle  of  the  freestream  boundary,  6max-  If  enough  points  were  placed  along  each  radial 
line,  then  the  solution  was  free  of  oscillations.  This  appeared  to  be  only  a  partial  fix  to  the 
problem,  as  this  method  was  very  dependent  on  how  many  nodes  were  placed  along  each 
radial  line.  The  method  that  gave  the  best  results  was  to  evaluate  the  continuity  equation 
at  node  k  for  the  interior  node  points,  and  at  the  midpoint  of  node  k  and  node  ks  for  the 
freestream  boundary  condition.  A  more  detailed  explanation  of  this  method  is  outlined  in 
Chapter  2.7. 


Wall 

Freestream 

u  =  0 

t;  =  0 

u>  =  0 
^=0 

p  =  1 

U  =  Uf, 

V=  Vf, 

W  =  Wf, 

e  =  ef. 

Table  4.1  Initial  boundary  conditions  for  viscous  subsonic  flow 


Wall 

Freestream 

1;=° 

V  =  0 

u;  =  0 
1^=0 

p=  1 

U  =  Uf, 

V  =  Vf, 

W  =  Wf, 

e  =  ef. 

Table  4.2  Initial  boundary  conditions  for  inviscid  subsonic  flow 


4-1  Supersonic  Results 

These  results  were  obtained  as  part  of  a  validation  procedure  for  the  analysis.  This 
investigation  provided  assurance  that  the  Jacobian  elements  were  being  placed  correctly 
in  the  Jacobian  matrix,  and  the  basic  logic  of  the  computer  code  was  correct.  Figure  4.1, 
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from  (2),  shows  the  pertinent  parameters  for  this  problem.  All  of  the  results  in  this  section 
were  obtained  for  a  =  0. 

For  supersonic,  inviscid,  axisymmetric  flow,  the  governing  equations  (2.15  -  2.19) 
reduce  to: 


2pu  sind  -|-  (pv  sin^)«  =  0, 

(4.1) 

p  [wiij  -  v^]  =  0, 

(4.2) 

P  [(y),  H  ^ 

(4.3) 

jn  =  0, 

(4.4) 

p[ve»\  -1-  (7  -  \)ep[2u  -1-  ucotd  -|-  ««]  =  0. 

(4.5) 

Compared  to  the  governing  equations  in  Chapter  2.3,  some  of  the  terms  in  the  equa¬ 
tions  are  cast  in  a  more  conservative  form  to  better  resolve  the  shock  position.  Several 
runs  were  made  at  varying  Mach  numbers  and  the  resulting  shock  positions  are  compared 
to  published  data  (28)  in  Table  4,1.  These  runs  were  made  with  six  radial  lines  around  the 
cone  and  either  201  or  401  node  points  equally  spaced  along  each  radial  line.  A  cone  angle, 
of  five  degrees  was  selected;  the  outer  radius  of  the  domain  was  specified  to  be  50  degrees 
ifimax  =  50).  For  all  cases  examined,  the  shock  angles  agree  very  well,  especially  when  401 
nodes  are  used.  The  shock  angle  is  determined  by  examining  the  non-dimensional  density 
between  the  cone  and  the  freestream  boundary.  As  shown  in  Figure  4.2,  the  density  varies 
from  a  value  greater  than  one  near  the  cone,  to  a  value  that  is  nearly  one  at  the  shock. 
The  shock  angle  is  determined  to  be  located  at  the  first  node  (starting  from  the  cone) 
where  the  density  becomes  one  (equal  to  the  freestream  non-dimensional  density).  For 
the  example  detailed  in  Figure  4.2,  the  Mach  number  is  1.321,  and  the  estimated  shock 
angle  is  49.375  deg.  With  401  nodes  used  in  the  radial  direction,  and  a  difference  in  and 
Omax  of  50  degrees,  the  angular  change  between  nodes  is  0.125  deg.  Therefore,  the  next 
computational  node  closer  to  the  cone  is  located  at  49.25  deg.  From  (22),  the  exact  shock 
angle  is  49.262  deg,  which  lies  in  between  the  two  computational  nodes  outlined  above.  If 
needed,  further  refinement  of  the  computational  domain  will  give  better  resolution. 
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Figures  4.3  and  4.4  show  the  computed  values  of  the  non-dimensional  velocity  com¬ 
ponents  u,  v  as  compared  to  the  exact  values,  computed  using  the  Taylor- MaccoU  equations 
and  tabulated  in  (22).  These  results  agree  extremely  well  over  the  entire  region  between 
the  shock  and  the  wall.  These  results  are  for  a  Mach  number  of  1.321,  chosen  to  match 
one  of  the  tables  in  (22). 

Subsonic  Results 

Several  point  solutions  were  attempted  in  the  development  of  the  solution  algorithm. 
The  effect  of  the  bandwidth  on  the  size  of  the  Jacobian  matrix  has  been  outlined  in 
Chapter  II.  As  the  bandwidth  increased,  the  computational  speed  decreased.  In  addition, 
increasing  the  bandwidth  also  increased  the  memory  requirements.  These  two  constraints 
imposed  a  limit  on  the  size  of  the  grid  that  could  be  used  in  the  analysis.  This  in  turn 
limits  the  parameters  used  in  a  point  solution.  Point  solutions  at  varying  angles  of  attack 
were  computed  for  the  flow  parameters  outlined  in  Table  4.6.  The  number  of  nodes  around 
the  cone  was  varied  between  30,  50,  and  81,  and  the  number  of  nodes  between  the  wall 
and  the  body  kept  at  81.  The  results,  presented  as  contour  plots  of  non-dimensional 
pressure,  are  shown  in  Figures  (4.5  -  4.14).  These  plots  are  correct  qualitatively,  although 
no  comparison  with  published  data  could  be  accomplished  due  to  the  low  Reynolds  number 
being  used  at  this  time.  Figure  4.15  shows  the  coefficient  of  pressure  at  the  cone  surface 
for  the  three  different  grids.  There  is  virtually  no  change  in  the  pressure  coefficient  as  the 
grid  is  changed,  indicating  that  at  least  for  this  case  the  boundary  layer  is  being  resolved 
correctly  and  the  solution  very  close  to  the  body  is  independent  of  the  grid.  The  values  of 
the  pressure  coefficient  as  a  function  of  <(>  are  consistent  with  published  data,  though  no 
direct  comparison  could  be  made  at  this  run  condition. 

At  first,  the  computer  used  for  this  project  was  only  capable  of  handling  a  30x81 
grid.  For  this  grid,  attempts  were  made  to  vary  the  amount  of  artificial  dissipation  added 
to  the  model.  It  was  expected  that  the  added  dissipation  would  help  smooth  the  pressure 
contours  away  from  the  body,  where  the  grid  coarseness  is  evident.  However,  as  more 
dissipation  was  added,  the  convergence  behavior  became  worse.  This  indicates  an  error  in 
the  implementation  of  the  artificial  dissipation,  which  has  not  been  discovered  at  this  time. 
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Solutions  at  higher  Reynolds  numbers  {Re  =  5000)  were  attempted,  but  the  convergence 
behavior  was  poor.  As  the  Reynolds  number  is  increased,  the  boundary  layer  becomes 
thinner.  This  requires  greater  resolution  near  the  wall,  which  means  there  will  be  less 
resolution  away  from  the  wall,  assuming  the  same  number  of  nodes  are  being  used.  Because 
of  the  constraint  on  the  number  of  nodes  due  to  the  memory  requirements,  increasing  the 
Reynolds  number  makes  it  harder  and  harder  for  the  solution  to  converge. 

At  this  point  a  different  computer  platform  became  available.  This  new  platform 
had  a  significant  increase  in  the  amount  of  available  memory,  as  well  as  computational 
speed.  The  code  was  implemented  on  this  machine,  and  grids  of  50x81  and  60x81  were 
now  possible.  Several  point  solutions  were  calculated,  using  the  parameters  outlined  in 
Table  4.6.  As  can  be  seen  from  the  pressure  contour  plots,  the  additional  nodes  around  the 
body  helped  to  smooth  the  contours  near  the  body,  but  there  is  stiU  evidence  of  spurious 
results  away  from  the  body.  At  these  grids,  a  point  solution  required  several  hours  to 
complete.  Therefore,  the  msun  effort  of  the  project  shifted  to  finding  ways  to  speed  up  the 
code.  With  the  implementation  of  a  block  structure  gauss  elimination  routine,  the  speed 
of  the  code  was  eventually  increased  by  a  factor  of  four  or  five.  EflTorts  were  also  made 
in  implementing  methods  that  would  decrease  the  memory  requirements  of  the  program, 
allowing  for  larger  grids  being  used.  These  modifications  were  only  recently  made,  and  the 
memory  savings  they  provide  could  not  be  taken  advantage  of  for  this  study. 

Addition  of  the  continuation  method  to  the  algorithm  was  accomplished.  With  this 
method,  a  point  solution  is  computed  at  zero  angle  of  attack  to  determine  the  sign  of  the 
determinate  of  the  Jacobian  matrix.  Then  continuation  occurs  in  an  attempt  to  locate 
the  bifurcation  point,  if  one  exists.  This  point  is  located  when  the  sign  of  the  determinate 
changes.  Once  the  bifurcation  point  is  located,  the  solution  method  is  perturbed  by  a 
small  value  so  the  nontrivial  solution  path  can  be  explored  (see  Figure  1.2a).  Only  manual 
continuation  was  attempted  in  this  study.  No  evidence  of  a  bifurcation  point  was  discov¬ 
ered,  but  without  a  grid  sensitivity  study,  there  is  no  way  to  determine  if  this  or  any  other 
result  is  accurate.  Also,  the  Reynolds  number  used  in  this  study  {Re  =  500),  may  be  to 
small  for  this  type  of  analysis. 
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It  should  be  noted  that  while  all  of  the  original  goals  of  this  project  were  not  com¬ 
pleted,  a  significant  amount  of  work  was  accomplished.  All  of  the  tools  required  for  the 
completion  of  the  analysis  have  been  developed.  The  computer  program  has  been  debugged 
and  validated.  The  continuation  method  was  implemented,  and  several  efforts  were  made 
to  increase  the  speed  of  the  code,  as  well  as  decrease  the  memory  requirements. 


4-6 


Table  4.3  Boundary  conditions  for  subsonic,  inviscid  flow 


WaU 

Freestream  Boundary 

continuity 

+  u’  -H  =  1 

p=  1 

W  =  Wft 

e  =  ef. 

Table  4.4  Boundary  conditions  for  subsonic,  viscous  flow 


Figure  4.1  Supersonic  Flow  over  Cone  at  a  =  0 
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Non-dimensional  density 


Figure  4.2  Non-dimensional  density 


Table  4.6  Run  Conditions  for  Point  Solutions 
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u  velocity 


Figure  4.3  Non-dimensional  u  velocity  component 
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V  velocity 


Figure  4.4  Non-dimensional  v  velocity  component 


Level  p 

D  7.98956 

C  7.98593 

B  7.98229 

A  7.97865 

9  7.97502 

8  7.97138 

7  7.96774 

6  7.9641 

5  7.96047 

4  7.95683 

3  7.95319 

2  7.94956 

1  7.94228 
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4-14 


Figure  4.10  50x81,  Re  =  500,  M  =  .3,a  =  8  deg 
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50x81  aoa  s  1 1  deg 


Level  p 

F  7.99803 

E  7.99305 

0  7.96804 

C  7.98302 

B  7.97801 

A  7.97299 

9  7.96798 

8  7.96296 

7  7.95794 

6  7.95293 

5  7.94792 

4  7.9429 

3  7.93789 

2  7.93287 

1  7.92786 


Figure  4.11  50x81,  Re  =  500,  M  =  .3,a  =  11  deg 
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50x81 
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5.0 

2.5 

0.0 

2.5 
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7.5 

-5  0 

Figure  4.12  50x81,  Re  =  500,  M 
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Level  p 
F  8.01858 

E  8.01103 

D  8.00348 

C  7.90593 

B  7.98838 

A  7.98083 

9  7.97328 

8  7.96573 

7  7.95818 

6  7.95063 

5  7.94308 

4  7.93553 

3  7.92798 

2  7.92043 

1  7.91288 


5 


10 


3,  a  =  15  deg 


60x81  aoa»11deg 


Lavoi  p 
F  8.01865 


Figure  4.14  60x81,  Re  =  500,  M  =  .3,a  =  15  deg 
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Figure  4.15 
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V.  Conclusions  and  Recommendations 


The  compressible,  laminar,  viscous  Navier-Stokes  equations  were  solved  for  subsonic 
flow  over  a  slender  circular  cone  at  angle  of  attack.  A  supersonic,  inviscid  study  was 
performed  at  zero  angle  of  attack  as  a  validation  of  the  solution  implementation.  Ex¬ 
cellent  result  were  obtained  for  the  supersonic  cases,  where  the  shock  angle  and  velocity 
components  were  calculated  and  compared  to  tabulated  results. 

Several  point  solutions  were  calculated  for  varying  angles  of  attack,  and  for  different 
grid  sizes.  A  continuation  method  was  implemented,  which  is  used  to  locate  the  bifurcation 
point,  if  one  exists.  For  this  investigation,  only  manual  continuation  was  used. 

The  most  obvious  problem  that  needs  to  be  analyzed  is  the  bandwidth  influence  on 
computational  speed  and  memory  requirements.  With  the  row-by-row  numbering  scheme, 
the  driving  requirement  is  the  use  of  second-order-accurate  approximations  for  the  bound¬ 
ary  derivatives  at  the  wall.  First-order-accurate  approximations  could  be  used,  but  the 
loss  of  accuracy  makes  this  undesirable.  An  alternative  solution  is  to  modify  the  computer 
code  so  that  the  analytical  Jacobian  elements  are  stored  in  single  precision.  When  any 
computations  are  made  with  these  elements,  double  precision  would  be  used.  This  effec¬ 
tively  cuts  the  memory  requirements  of  the  computer  code  in  half.  There  is  some  loss  of 
accuracy  and  convergence  rate,  but  these  disadvantages  should  be  tolerable. 

Recent  analysis  has  dealt  with  speeding  up  the  code  by  modifying  the  structure  of 
the  matrix  so  that  the  system  of  equations  could  be  solved  using  a  blocked  method.  This 
has  been  implemented  and  the  result  is  a  significant  speed  up  of  the  code.  Unfortunately, 
this  does  nothing  for  the  memory  issue. 

Much  more  effort  is  required  to  complete  all  the  analysis  necessary  to  determine 
if  asymmetric  vortices  are  a  real  solution  of  the  Navier-Stokes  equations.  An  extensive 
study  of  the  grid,  and  its  effect  on  the  solution  needs  to  be  performed.  Areas  of  interest 
include,  but  are  not  limited  to,  grid  refinement  in  both  the  <f>  and  9  directions,  choosing 
an  amount  of  artificial  dissipation  that  allows  for  a  smooth  solution  without  destroying 
the  accuracy  of  the  results  near  the  body,  examining  the  effect  of  the  boundary  conditions 
on  the  results,  and  examining  the  boundary  layer  resolution  by  investigating  solutions 


at  different  Reynolds  numbers.  Analysis  of  this  type  will  help  to  characterize  the  length 
scales  in  the  problem,  and  deternnine  whether  a  particular  grid  will  resolve  all  of  these 
length  scales.  Once  a  set  of  parameters  is  determined  that  will  result  in  a  grid  independent 
solution,  continuation  can  be  used,  either  manual  or  implemented  in  the  program,  to  vary 
the  angle  of  attack  in  an  attempt  to  locate  a  bifurcation  point.  If  this  step  is  reached,  the 
solution  method  implemented  in  this  computer  model  can  be  used  to  describe  all  branches 
in  the  solution  space,  whether  stable  or  unstable. 
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Appendix  A.  Derivation  of  Governing  Equations 


The  derivation  of  the  non-dimensional  governing  equations  begins  with  the  dimen¬ 
sional  Navier-Stokes  equations  of  motion  for  a  compressible  fluid  in  spherical  coordinates, 
as  found  in  reference  (17).  The  five  governing  equations  consist  of  conservation  of  mass, 
momentum  (three  equations),  and  internal  energy: 
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where  qr  is  the  radiation  heat  flux  vector  and  is  assumed  to  be  zero.  Q  represents  internal 
heat  generation  and  is  also  assumed  to  be  zero.  The  body  forces  Fr,Ft,F^  are  assumed 
negligible  and  so  are  discarded.  $  is  the  viscous  dissipation  function: 
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The  material  derivative  in  spherical  coordinates  is  defined  as 
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where  V  is  the  velocity  vector.  Two  of  the  major  assumptions  in  this  study  are  steady- 
state  flow  and  conical  similarity.  Steady  state  implies  |f  =  0  and  conical  similarity  implies 
=  0,  where  /  is  one  of  the  dependent  variables.  The  pressure,  P,  is  written  in  terms 
of  the  density,  p,  and  the  internal  energy,  e,  assuming  a  thermally  and  calorically  perfect 
gas,  by 
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Elimination  of  T  from  (A.9)  yields 
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The  temperature,  T,  is  written  in  terms  of  e  using  equation  (A. 10).  Stokes’  hypothesis 
is  also  used,  A  =  The  equations,  after  all  terms  are  expanded  and  the  above  re¬ 

lationships  implemented,  are  then  non-dimensionalized.  The  dependent  and  independent 
variables  are  non-dimensionalized  as  follows: 
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The  coordinate  directions  0  and  </>  are  already  non-dimensional  angles,  so  they  are 
not  affected  in  this  process.  The  non-dimensional  viscosity  is  replaced  with  the  linear 
relationship  p,  =  CjC  +  C2,  developed  in  Chapter  2,  providing  the  final  form  of  the  non- 
dimensional  equations  (after  dropping  the  (*)  for  convenience): 

Continuity  equation 
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9-momentum  equation 
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</>-momentum  equation 
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Energy  equation 
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(A.15) 


(A.16) 


(A.17) 
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To  check  to  accuracy  of  the  preceeding  derivation,  the  continuity  and  momentum 
equations  were  re-derived  using  the  non-dimensional  Navier-Stokes  equations  found  in 
Reference  (31).  The  continuity  equation  is 


The  momentum  equations  for  the  ii,  *21  a-^d  *3  coordinate  directions  are,  respectively, 
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The  components  of  the  tensor  T,  Tij ,  are  given  by 


(A.21) 


Tij  =  pUiUj  -1-  P6ij  -  (A.22) 

where 
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(A.27) 


(A.28) 


/J  1  ,  Uidha  .  Ujd/ia^ 

(J./V)„=-(^— +  ^_  +  -  — j  (A.27) 

and  Sij  is  the  Kronecker  delta.  For  a  spherical  coordinate  system, 

hi  =  1  Xi  =  r  Ui  =  u 

h2  =  r  X2  =  6  U2  =  V 

hs  =  rsind  13  =  <f>  U3  =  w. 

When  equations  (A.  19  -  A.22)  are  expanded  and  the  continuity  equation  (A. 19) 

subtracted  out  of  the  momentum  equations,  the  results  are  the  same  as  the  previous 
derivation  for  (A. 15  -  A.17).  The  energy  equation  was  not  checked  through  an  alternate 
derivation. 
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Appendix  B.  Derivation  of  Analytical  Jacobian  Elements 

As  an  example  of  how  all  the  equations  are  manipulated,  the  continuity  equation  will 
be  transformed  from  (^,  9)  physical  space  to  (^,  q)  computational  space  and  the  analytical 
Jacobian  elements  derived.  The  other  four  equations  are  then  handled  in  the  same  manner 
to  determine  their  respective  Jacobian  elements. 

With  the  subscripts  9  and  4>  representing  differentiation,  the  non-dimensional  conti¬ 
nuity  equation  in  spherical  coordinates,  as  developed  in  Appendix  A,  becomes 

2pu  sind  -f-  pv  CO&9  -1-  vpt  sind  -i-  pvt  sin^  -t-  pw^  -I-  wp^  =  0. 

The  derivatives,  transformed  using  the  relationships  outlined  in  Chapter  2,  are 

P<t>  —  i>t>Pi  +  'n<^Pni  P»  =  ^0Pi  +  V»Pn 


The  ^  and  rj  derivatives  are  now  expressed  in  terms  of  the  nine-'  ^mt  stencil,  shown 
in  Figure  2.4,  using  central  differences  about  node  k.  The  central  differences  needed  for 
the  continuity  equation  are 
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Approximations  (B.l  -  B-3)  are  second-order-accurate  in  the  node  spacing.  Using  a  square 
computational  cell,  =  At/  =  1,  the  continuity  equation  becomes 
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Figure  B.l  Nine-point  Stencil  for  interior  nodes 
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The  analytical  Jacobian  elements  are  now  determined  by  taking  partial  derivatives 
of  (B.4)  with  respect  to  each  of  the  dicrete  flow  variables  in  the  nine-point  stencil,  shown 
in  Figure  B.l.  The  Jacobian  elements  for  the  continuity  equation  are  found  to  be 
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